1 analysis set up

1.1 load libraries and source used functions

1.2 read in gene count table

Trimmed and filtered reads were aligned using STAR (v2.7.0f). Read counting was performed using featureCounts (v1.6.2) with parameters “-s 1 -Q 255”. Gene annotation was generated by merging GENCODE (vM21), 3’UTR annotations in RefSeq (mm10, UCSC), and de novo 3’UTR sites discovered from the data.

Count table is read into R. A sample condition table is generated by parsing filenames.

# gene count file
counts_file <- "datasets/pit_utrseq_gencode_vM21_refSeq201906.txt"

# Obtain gene count table and sample conditions
temp <- get_count_table(counts_file, fixname = T)
count_table <- temp[[1]]
sampleCondition <- temp[[2]]
gene_info <- temp[[3]]
rm(temp)

colnames(count_table) <- gsub("p", "P", colnames(count_table))

# Get all chromosome Y genes and 3 chrX genes that are consistently female-biased
chrXYgenes <- as.character(gene_info[grepl("chrY", gene_info[,2]),1])
chrXYgenes <- c(chrXYgenes, "ENSMUSG00000086503.3", "ENSMUSG00000086503.4", "ENSMUSG00000099312.1")
chrMgenes <- as.character(row.names(gene_info[grepl("chrM", gene_info[,2]),]))

2 Normalization using RUVseq

We used RUVseq to remove unwanted variables with empirically identified negative control genes.

# constructc edger expression set to obtain cpm values
y <- DGEList(counts=count_table, genes=row.names(count_table), group=sampleCondition$group)
y <- calcNormFactors(y)
count_cpm <- cpm(y)

# get spike-in gene names
spikes <- rownames(count_table)[grep("^ERCC", row.names(count_table))]

# filter for only genes expressed at cpm > 2 in at least 10 samples. Also removed mitochondrial genes and spike-ins
count_table_fil <- count_table[which(rowSums(count_cpm > 2) >= 10), ]
count_table_fil <- count_table_fil[!row.names(count_table_fil) %in% chrMgenes,]
count_table_fil <- count_table_fil[!row.names(count_table_fil) %in% spikes,]

# construct edger expression set again with filtered data
y_fil <- DGEList(counts=count_table_fil, genes=row.names(count_table_fil), group=sampleCondition$age)
y_fil <- calcNormFactors(y_fil)

# get log CPM counts
logCPM <- cpm(y_fil, log=TRUE, prior.count=1)


# Correct batch effect/unwanted variables using RUV seq 
 
row.names(sampleCondition) <- colnames(count_table_fil)
set <- newSeqExpressionSet(as.matrix(count_table_fil),
                           phenoData = sampleCondition)
# upper quartile normalzation
set <- betweenLaneNormalization(set, which="upper")


# identify empirical control genes by performing an DE analysis across all conditions
design <- model.matrix(~group, data=pData(set))
y_fil <- DGEList(counts=count_table_fil, genes=row.names(count_table_fil), group=sampleCondition$group)
y_fil <- calcNormFactors(y_fil)
y_fil <- estimateDisp(y_fil, robust = TRUE, design = design)

fit <- glmQLFit(y_fil, robust = T, design)
lrt <- glmQLFTest(fit, coef=2:10)

top <- topTags(lrt, n=nrow(count_table_fil))$table
top$genename <- genename[top$genes]

# identified any genes with an FDR > 0.1 as negative controls
empirical <- subset(top, FDR > 0.1)$genes

# RUVseq normliazation
set1 <- RUVg(set, empirical, k=2)

# obtain normalized counts
logCPMc <- log2(normCounts(set1)+1)

3 Overview

3.1 Comparing quantification of 3’UTRseq and qPCR (Figure 1C, Supp Figure 1C)

qPCR of 183 puberty-related and control genes have been performed previously using the same samples. Here, we compared the expression quantification using 3’UTRseq (log2 normalized read counts) and that using qPCR (delta Ct)

# get qPRC results
all_data_normed_info_noC <- readRDS("datasets/fluidigm_all_data_normed_info_noC.rds")
# fix gene names
all_data_normed_info_noC$Primer[all_data_normed_info_noC$Primer == "Rab7l1"] <- "Rab29"
fluidigm_genes <- unique(all_data_normed_info_noC$Primer)
all_genes <- fluidigm_genes

# reformat pcr results table
PCR_puberty <- subset(all_data_normed_info_noC, tissue == "P") %>%
  mutate(Sample = gsub("_", "", Sample)) %>%
  group_by(Sample, Primer, age, gender) %>%
  summarise(Ct= mean(Ct)) %>%
  as.data.frame()

# compare pcr results and utrseq results (before and after normalization, only post-normalization data shown)
exp_compare_cor_values <- compare_qPCR(logCPM, logCPMc, PCR_puberty, outname, method = "spearman", ifgenplot = F)

compare_sum <- compare_qPCR_summaryplot(exp_compare_cor_values$cor_values)

exp_compare_cor_values$plots_after_correction$Pd37M4
Figure S1C: summary of UTRseq vs. qPCR quantification comparison

Figure 1: Figure S1C: summary of UTRseq vs. qPCR quantification comparison 

compare_sum[[2]]
Figure 1C: example scatter plot comparing UTRseq and qPCR

Figure 2: Figure 1C: example scatter plot comparing UTRseq and qPCR

Pearson correlation coefficient between normalized UTRseq counts and qPCR delta Ct ranges from 0.7 to 0.78.

3.2 Sample global overview (Figure 1D and Supp Figure 2A )

PCA plot and correlation heatmap using RUVseq normalized counts (log2 transformed).

# get PCA plots
logCPMc_PCAplot <- gen_pca_plot(logCPMc, npc=20, sample_label=F, label =F)
logCPMc_PCAplot[[1]]
Figure 1D: PCA plot of all samples

Figure 3: Figure 1D: PCA plot of all samples 

hmcol <- colorRampPalette(c("orange", "white", "darkmagenta"))(250)
plotHM(logCPMc, cols = hmcol, new_col_names = sapply(strsplit(colnames(logCPMc), "_"), "[[", 2), show_col_names = F, clust_method = "complete", anno = sampleCondition)
Figure S2A: correlation heatmap of samples

Figure 4: Figure S2A: correlation heatmap of samples 

# calculate mean person correlation between biological reps
sample_cors <- cor(logCPMc)
sample_cons <- unique(sampleCondition$group)
x <- sample_cons[1]
sample_cors_means <- lapply(sample_cons, function(x) {
  cors <- sample_cors[grepl(x, row.names(sample_cors)), grepl(x, colnames(sample_cors))]
  cors[cors==1] <- NA
  return(list(mean(cors, na.rm=T), cors))}
  )
mean_cor <- round(sapply(sample_cors_means,"[[", 1),2)

Average pearson correlation between biological replicates (samples from the same biological condition) ranges from 0.95 to 0.97

4 DE analysis

Pairwise differential analysis was performed using EdgeR. Genes with a FDR < 0.05 and an absolute fold change> 1.5 are considered significant.

Comparisions are performed between male and female samples at each age

4.1 DE comparison between sex

all_detected_genes <- genename[row.names(count_table_fil)]
y_fil <- DGEList(counts=count_table_fil, genes=row.names(count_table_fil), group=sampleCondition$group)
y_fil <- calcNormFactors(y_fil)

design <- model.matrix(~0+group+W_1+W_2, pData(set1))

y_fil <- estimateDisp(y_fil, robust = TRUE, design = design)
fit <- glmQLFit(y_fil, design=design, robust = TRUE)


# compare between sexes at each age.
pair_contrasts <- makeContrasts(d12_sex = groupd12F - groupd12M,
                                  d22_sex = groupd22F - groupd22M,
                                  d27_sex = groupd27F - groupd27M,
                                  d32_sex = groupd32F - groupd32M,
                                  d37_sex = groupd37F - groupd37M,
                                  levels = design)
  
  
  get_test_result <- function(contrasts=pair_contrasts, name){
    qlf <- glmQLFTest(fit, contrast = contrasts[, name])
    topTg <- topTags(qlf, n=nrow(y_fil$counts))
    
    de <- as.data.frame(topTg[[1]])
    de$genename <- genename[de$genes]
    return(de)
  }
  
de_result_list <- lapply(colnames(pair_contrasts), function(x) get_test_result(pair_contrasts, x))
names(de_result_list) <- colnames(pair_contrasts)
de_sig_list <- lapply(de_result_list, function(x) subset(x, FDR < 0.05 & abs(logFC) > log2(1.5)))
names(de_sig_list) <- colnames(pair_contrasts)

templist <- flatten(lapply(de_sig_list, function(x) list(subset(x, logFC>0), subset(x, logFC<0))))

names(templist) <- flatten(lapply(names(de_sig_list), function(x){
  sapply(c("up", "down"), function(y) paste(x, y, sep = "_"))
}))

sig_genes_names_list_direction <- lapply(templist, "[[", "genename")
sig_genes_names_list <- lapply(de_sig_list, "[[", "genename")

4.1.1 Overlap between sex-biased genes (Figure 2A)

sex_sig_genes_df <- get_sig_genes_df_from_list(de_sig_list)

colnames(sex_sig_genes_df) <- gsub("_sex_up", "_F", colnames(sex_sig_genes_df))
colnames(sex_sig_genes_df) <- gsub("_sex_down", "_M", colnames(sex_sig_genes_df))

upset(sex_sig_genes_df, sets=rev(c("d12_F", "d22_F", "d27_F", "d32_F","d37_F", "d12_M", "d22_M", "d27_M","d32_M", "d37_M")), keep.order = T,
      order.by = "freq", 
      queries = list(list(query = intersects, params = list("d27_F", "d32_F", "d37_F"), color="tomato", active=T),
                     list(query = intersects, params = list("d27_M","d32_M", "d37_M"), color="steelblue", active=T)))
Figure 4A: overlap between sex-biased genes

Figure 5: Figure 4A: overlap between sex-biased genes 

# get sex biased genes in 2 ages
temp <- table(unlist(sig_genes_names_list_direction[c("d27_sex_up", "d32_sex_up","d37_sex_up")]))
f_biased_genes_used <- names(temp[temp !=1])
f_biased_genes_3 <- names(temp[temp ==3])

temp <- table(unlist(sig_genes_names_list_direction[c("d27_sex_down", "d32_sex_down","d37_sex_down")]))
m_biased_genes_used <- names(temp[temp !=1])
m_biased_genes_3 <- names(temp[temp ==3])

# f_biased_genes_used and m_biased_genes_used are used as input for Lisa

4.1.2 selecting example sex-biased genes (Figure 2C-D)

sig_df <- melt(sig_genes_names_list) %>% 
  separate(L1,into = c("age",NA))

# plot selected d12-d22 sex biased genes
ave_expr <- rowMeans(logCPMc)
names(ave_expr) <- genename[names(ave_expr)]

# figure2 genes are selected based on 1) d12 sex biased genes , 2) gene function mentioned in the text and 3) expression levels 
figure2_genes1 <- c("Chrna4","Kcna4","Th","Drd4", "Asb4", "Fshb","Lhb","Serpine2", "Nupr1","Gpr101")

factor_plot(logCPMc, figure2_genes1, ncols = 5, sig_df = sig_df)

# get top m-/f-biased genes that are consistent from PD27
top_fbiased <- sort(ave_expr[names(ave_expr) %in% f_biased_genes_3], decreasing = T)[1:5]
top_mbiased <- sort(ave_expr[names(ave_expr) %in% m_biased_genes_3], decreasing = T)[1:5]

factor_plot(logCPMc, c(names(top_fbiased), names(top_mbiased)), ncols = 5, sig_df = sig_df)

4.1.3 Pathway enrichment of sex-biased genes

enrich_list_direction <- lapply(sig_genes_names_list_direction, function(x) gprofiler(x, organism = "mmusculus", custom_bg = all_detected_genes, hier_filtering = "moderate", src_filter = c("GO", "KEGG","REAC","CORUM","HP","HPA","OMIM"),min_set_size = 5, max_set_size = 5000,min_isect_size = 3))

male_biased_pathways <- enrich_list_direction[c("d27_sex_down", "d32_sex_down", "d37_sex_down")]
female_biased_pathways <- enrich_list_direction[c("d27_sex_up", "d32_sex_up", "d37_sex_up")]

maleenrichmentmap <- enrichmap_like_pie(male_biased_pathways, layout_used = "kk", radius_scale = 0.05)
femaleenrichmentmap <- enrichmap_like_pie(female_biased_pathways, layout_used = "kk", radius_scale = 0.07)

4.1.3.1 pathway enrichment of male-based genes (Figure 2F)

maleenrichmentmap[[2]]
Figure 2F: enrichment map-like representation of pathway enrichment of male genes at PD27, 32, and 37

Figure 6: Figure 2F: enrichment map-like representation of pathway enrichment of male genes at PD27, 32, and 37 

maleenrichmentmap[[3]]
Figure 2F: enrichment map-like representation of pathway enrichment of male genes at PD27, 32, and 37

Figure 7: Figure 2F: enrichment map-like representation of pathway enrichment of male genes at PD27, 32, and 37

4.1.3.2 pathway enrichment of female-based genes (Figure 2G)

femaleenrichmentmap[[2]]
Figure 2G: enrichment map-like representation of pathway enrichment of female genes at PD27, 32, and 37

Figure 8: Figure 2G: enrichment map-like representation of pathway enrichment of female genes at PD27, 32, and 37 

femaleenrichmentmap[[3]]
Figure 2G: enrichment map-like representation of pathway enrichment of female genes at PD27, 32, and 37

Figure 9: Figure 2G: enrichment map-like representation of pathway enrichment of female genes at PD27, 32, and 37

4.1.4 Testing for sex-age interaction

# sex age interaction
sex_age_contrasts <- makeContrasts(d12_d22 = (groupd22F - groupd12F) - (groupd22M - groupd12M),
                                   d22_d27 = (groupd27F - groupd22F) - (groupd27M - groupd22M),
                                   d27_d32 = (groupd32F - groupd27F) - (groupd32M - groupd27M),
                                   d37_d32 = (groupd37F - groupd32F) - (groupd37M - groupd32M),
                                   d32_d22 = (groupd32F - groupd22F) - (groupd32M - groupd22M),
                                   d37_d27 = (groupd37F - groupd27F) - (groupd37M - groupd27M),
                                   d37_d22 = (groupd37F - groupd22F) - (groupd37M - groupd22M),
                                   levels = design)


de_sex_age_result_list <- lapply(colnames(sex_age_contrasts), function(x) get_test_result(sex_age_contrasts, x))
names(de_sex_age_result_list) <- colnames(sex_age_contrasts)
de_sex_age_sig_list <- lapply(de_sex_age_result_list, function(x) subset(x, FDR < 0.05 & abs(logFC) > log2(1.5)))

sig_sex_age_genes_names_list <- lapply(de_sex_age_sig_list, function(x) x$genename)
#sig_genes_puberty <- lapply(sig_genes_names_list, function(x) x[x %in% all_puberty_genes])

d12_d22_sex_diff <- de_sex_age_sig_list$d12_d22$genename
#factor_plot(logCPMc, sort(d12_d22_sex_diff), ncols=5, name=F)

4.1.5 Heatmap showing the expression of selected sex-biased genes (Figure S3A)

The genes shown in this heatmap includes:

  • all sex-biased genes at PD12 and PD22
  • genes with PD12-PD22 age-sex interaction effect

that are not already shown in Figure 2

# get sex biased genes in 2/3 ages

sex_biased_genes_to_show <- unique(c(unlist(sig_genes_names_list[1:2]),d12_d22_sex_diff))

# exclude genes that are shown in figure 2 already
sex_biased_genes_to_show <- sex_biased_genes_to_show[!sex_biased_genes_to_show %in% c(figure2_genes1, names(top_fbiased), names(top_mbiased))]

#sex_biased_genes_to_show <- c(sex_chr_genes, "Kcna4","Th", d22_rest)

sex_genes_sum <- table(unlist(sig_genes_names_list[1:5]))

temp_df <- data.frame(value = sex_biased_genes_to_show, stringsAsFactors = F)

sex_gene_df <- left_join(temp_df, melt(sig_genes_names_list_direction)) %>% 
# filter(value %in% sex_biased_genes_to_show) %>% 
  separate(L1, "_", into = c("age", NA, "dir")) %>% 
  spread(age, dir) %>% 
  mutate(n = sex_genes_sum[as.character(value)]) %>% 
  dplyr::select(value, d12, d22, d27, d32, d37) %>% 
  arrange(d12, d22, d27, d32, d37) %>% 
  column_to_rownames("value")

#sex_gene_df$d12_d22 <- ifelse(row.names(sex_gene_df) %in% d12_d22_sex_diff, "yes","no")

anno_cols_sex <-rep(list(c("up" = "tomato","down" = "steelblue")), 5)
names(anno_cols_sex) <- c("d12","d22","d27","d32","d37")

#anno_cols_inter <- list("d12d22" = c("yes" = "darkmagenta"))

p <- plot_gene_hm(row.names(sex_gene_df), 
                  logCPMc, 
                  clust_cols = F, 
                  clust_rows = F,
                  row_cut = 1, 
                  anno_colors = c(anno_cols, anno_cols_sex), 
                  anno_rows = sex_gene_df[, c("d12","d22","d27","d32","d37")],
                  colnames = F,
                 # fontsize_row = 8,
                  colors = hmcol_yp, 
                  sampleCon = sampleCondition)


plot(p$gtable)
Figure S3A: heatmap of selected sex-biased genes

Figure 10: Figure S3A: heatmap of selected sex-biased genes

4.1.6 Expression heatmap of sex-biased TFs (Figure 3B)

These genes are TFs that are sex-biased at 2 out of 3 ages (PD27, 32, and 37) and are not shown in Figure 3.

# plots sex-biased tfs
sex_tfs <- c(m_biased_genes_used[m_biased_genes_used %in% alltfs], f_biased_genes_used[f_biased_genes_used %in% alltfs])

temp_df <- data.frame(value = sort(sex_tfs), stringsAsFactors = F)

sex_gene_df <- left_join(temp_df, melt(sig_genes_names_list_direction)) %>% 
  separate(L1, "_", into = c("age", NA, "dir")) %>% 
  filter(!is.na(age)) %>% 
 # spread(age, dir) %>% 
  pivot_wider(names_from = "age", values_from = "dir") %>% 
  mutate(sex = ifelse(d27=="down"|d32=="down"|d37=="down","down","up")) %>% 
  group_by(sex) %>% 
  arrange(d27, d32, d37,.by_group = T) %>% 
  column_to_rownames("value")


anno_cols_sex <-rep(list(c("up" = "tomato","down" = "steelblue")), 5)
names(anno_cols_sex) <- c("d12","d22","d27","d32","d37")

p_hm_tf <- plot_gene_hm(row.names(sex_gene_df),
                  logCPMc, 
                  clust_cols = F, 
                  clust_rows = F,
                  row_cut = 1, 
                  anno_colors = c(anno_cols, anno_cols_sex), 
                  anno_rows = sex_gene_df[, c("d27","d32","d37")],
                  colnames = F,
                 # fontsize_row = 8,
                  colors = hmcol_yp, 
                  sampleCon = sampleCondition)

plot(p_hm_tf$gtable)
Figure 3B: Expression of sex-biased TFs

Figure 11: Figure 3B: Expression of sex-biased TFs

4.2 testing if sex-biased genes are enriched forpituitary disease and pubertal timing associated genes

run_hyper <- function(query, geneset, total){
  q <- length(intersect(query, geneset)) #num overlap geneset
  m <- length(intersect(geneset, total)) # num of geneset
  n <- length(total[!total %in% geneset]) # background that are not geneset
  k <- length(intersect(query, total)) # num query genes
  
  fold <- round(q/(k*m/(m+n)),3)
  p <- phyper(q, m, n, k, lower.tail = F)
  return(c("num_ol"=q,"num_genes_in_geneset"=m,"num_background_noneGeneset"=n,"num_query"=k,"fold_enrichment"=fold,"pval"=p))
}

# read in all disease genes 
disease_genes <- read.table("datasets/pituitary_puberty_genes_combined_2020-09-28.txt", header = T, sep = "\t", stringsAsFactors = F, as.is = T)
all_disease_genes <- disease_genes$MGI.symbol

all_sex_biased_genes <- unique(unlist(sig_genes_names_list_direction[1:10]))

hyper_result <- run_hyper(all_sex_biased_genes, all_disease_genes, all_detected_genes)

hyper_result[c("fold_enrichment","pval")]
## fold_enrichment            pval 
##    1.801000e+00    2.667187e-05

5 Co-expression analysis

5.1 building co-expresison modules using CEMtools

# run CEMtools with default parameters
cem <- cemitool(as.data.frame(logCPMc), merge_similar = T, network_type = "signed")

# get module assignment
modules <- module_genes(cem)
modules$genename <- genename[modules$genes]
cemi_genes <- modules$genes

# transform into a list format
module_list <- lapply(split(modules, modules$modules), "[[", "genes")
module_list_genename <- lapply(split(modules, modules$modules), "[[", "genename")

# get average scaled expression for each gene for plotting in the next section
ave <- ave_count_table(logCPMc)

# get top hub genes
hubs <- CEMiTool::get_hubs(cem, n=20)
hubs_genename <- lapply(hubs, function(x) genename[names(x)])

5.2 co-expression module summary (Figure 4)

-Left:Heatmap showing the scaled normalized expression levels of each gene, clustered separately within modules.

-Middle: plot summarizing scaled normalized expression patterns of each gene in each module. Thick line represent the average profile of the module.

-Right: expression profile of the top hub genes of each module.

-Bottom: top 10 hub genes for each module

# scale normalized count data (log transformed) and reorder columns
plotdata <- t(scale(t(logCPMc)))
colnames(plotdata) <- sapply(strsplit(colnames(plotdata), "_"), "[[", 2)
plotdata <- plotdata[, c(colnames(plotdata)[grepl("M", colnames(plotdata))], colnames(plotdata)[grepl("F", colnames(plotdata))])]


get_clust_gene_order <- function(clust_num){
  # function to cluster genes within each cluster
  used_genes <- module_list[[clust_num]]
  used_data <- plotdata[row.names(plotdata) %in% used_genes, ]
  phm <- pheatmap(used_data, 
                  clust_cols = F,
                  clustering_distance_rows = "correlation",
                  clustering_method = "ward.D2",
                  silent = T)
  order <- used_genes[phm$tree_row$order]
  return(order)
}

# reorder genes based on clustering results within each cluster
reorder <- unlist(lapply(1:9, get_clust_gene_order))
cluster_length <- sapply(module_list, length)

# column annotation data
col_anno <- pData(set1)
row.names(col_anno) <-sapply(strsplit(row.names(col_anno), "_"), "[[", 2)

# set cutoffs for extreame values
all_values <- as.vector(as.matrix(plotdata))
custome_breaks <- unique(c(min(all_values), seq(quantile(all_values, 0.005), 0, length = 125), seq(0, quantile(all_values, 0.995), length=125), max(all_values)))

# generate heatmap
phm_modulegenes <- pheatmap(plotdata[reorder, ],
         cluster_rows = F, 
         cluster_cols = F, 
         show_rownames = F, 
         annotation_col = col_anno[, c("age", "sex")], 
         annotation_colors = anno_cols,
         color = hmcol, 
         breaks = custome_breaks,
         gaps_row = cumsum(cluster_length),
         silent = T)

# get a summarized plot of module gene expression
ave <- ave_count_table(logCPMc)

module_assignment <- modules$modules
names(module_assignment) <- modules$genes

pclusters <- plot_clusters(module_assignment, ave, "", ncols=1) + theme(axis.text.x = element_blank(), axis.text.y = element_text(size=14), strip.text = element_text(size=14))

# plot top module hub genes
top_hub_genes <- c("Tshb", "GH", "Prl", "Fn1", "Dlk1", "Ank1", "Mki67", "Col23a1", "Fshb")
tophub_genes <- factor_plot(logCPMc, top_hub_genes, ncols = 1)

hub_df <- as.data.frame(sapply(hubs_genename, function(x) paste(x[1:10], collapse = ",")))
colnames(hub_df) <- "Top10_hub_genes"

# put together
p <- as_ggplot(phm_modulegenes$gtable) + pclusters + tophub_genes +
  plot_layout(widths = c(4,1.5,1)) +
  theme_bw()

p + gridExtra::tableGrob(hub_df)
Figure 4:co-expression module expression and hub genes

Figure 12: Figure 4:co-expression module expression and hub genes

5.3 More example hub genes (Figure S4)

top5_hub_genes <- unname(unlist(lapply(hubs_genename, function(x) x[1:5])))

factor_plot(logCPMc, top5_hub_genes, ncols = 5)
Expression patterns of top 5 hub genes

Figure 13: Expression patterns of top 5 hub genes

5.4 Module correlation (Figure S5A)

module_ME <- mod_summary(cem, "eigengene")
module_ME_df <- melt(module_ME) %>%
  separate(variable, "_", into=c("id", "condition")) %>%
  mutate(age=substr(condition, 2,4), sex=substr(condition, 5,5), rep=substr(condition, 6,6)) 
row.names(module_ME) <- module_ME$modules
module_ME <- module_ME[, -1]

p <- pheatmap(cor(t(module_ME)), silent = T)
plot(p$gtable)
Figure S5A: correlation of module eigengenes

Figure 14: Figure S5A: correlation of module eigengenes

5.5 Module eigengene (Figure S5B)

ggplot(module_ME_df) +
  geom_rect(xmin = -Inf,xmax = Inf,
            ymin = -Inf,ymax = Inf,fill="gray90", alpha = 0.1) +
  geom_bar(aes(x=id, y=value, fill=sex), stat="identity") +
  facet_grid(modules~age, scales="free") +
  # theme(axis.text.x = element_text(angle=45, hjust=1)) +
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank()) +
  scale_fill_manual(values = c("tomato", "steelblue"))
Figure S5B: module eigengenes

Figure 15: Figure S5B: module eigengenes

5.6 Pathway enrichment analysis of module genes (Figure S5C)

module_gene_enrich <- lapply(module_list_genename, function(x) gprofiler(x, organism = "mmusculus", custom_bg = all_detected_genes, hier_filtering = "moderate", min_set_size = 5, max_set_size = 5000, min_isect_size = 3, src_filter = c("GO")))
p1 <- plot_cluster_enrichment(enrich_list = module_gene_enrich[1:4], "CEMi_modules_enrichment_", width=13, exclude_domain = "hp", save = F)
p2 <- plot_cluster_enrichment(enrich_list = module_gene_enrich[5:9], "CEMi_modules_enrichment_", width=13, exclude_domain = "hp", save = F)

p1+p2 +
  plot_layout(guides = "collect")
Figure S5C: pathway enrichment result of co-expression modules

Figure 16: Figure S5C: pathway enrichment result of co-expression modules 

save.image("mouse_pituitary_mRNA_analysis.RData")

6 SessionInfo

sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.3 LTS
## 
## Matrix products: default
## BLAS:   /opt/R/3.6.0/lib/R/lib/libRblas.so
## LAPACK: /opt/R/3.6.0/lib/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_CA.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_CA.UTF-8        LC_COLLATE=en_CA.UTF-8    
##  [5] LC_MONETARY=en_CA.UTF-8    LC_MESSAGES=en_CA.UTF-8   
##  [7] LC_PAPER=en_CA.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] pheatmap_1.0.12             CEMiTool_1.8.2             
##  [3] UpSetR_1.4.0                patchwork_1.0.0            
##  [5] scatterpie_0.1.4            ggraph_1.0.2               
##  [7] tidygraph_1.1.2             gProfileR_0.6.7            
##  [9] fgsea_1.10.0                Rcpp_1.0.1                 
## [11] dendextend_1.12.0           ggdendro_0.1-20            
## [13] cowplot_0.9.4               gridExtra_2.3              
## [15] reshape2_1.4.3              RUVSeq_1.18.0              
## [17] EDASeq_2.18.0               ShortRead_1.42.0           
## [19] GenomicAlignments_1.20.1    SummarizedExperiment_1.14.0
## [21] DelayedArray_0.10.0         matrixStats_0.54.0         
## [23] Rsamtools_2.0.0             Biostrings_2.52.0          
## [25] XVector_0.24.0              BiocParallel_1.18.0        
## [27] ggrepel_0.8.1               ggpubr_0.2.1               
## [29] magrittr_1.5                forcats_0.4.0              
## [31] stringr_1.4.0               dplyr_0.8.3                
## [33] purrr_0.3.3                 readr_1.3.1                
## [35] tidyr_1.0.0                 tibble_2.1.3               
## [37] ggplot2_3.2.1               tidyverse_1.2.1            
## [39] pcaMethods_1.76.0           Biobase_2.44.0             
## [41] gplots_3.0.1.1              RColorBrewer_1.1-2         
## [43] rtracklayer_1.44.4          GenomicRanges_1.36.1       
## [45] GenomeInfoDb_1.20.0         IRanges_2.18.1             
## [47] S4Vectors_0.22.0            BiocGenerics_0.30.0        
## [49] edgeR_3.26.5                limma_3.40.2               
## 
## loaded via a namespace (and not attached):
##   [1] ggthemes_4.2.0         R.methodsS3_1.7.1      coda_0.19-2           
##   [4] ggpmisc_0.3.1          acepack_1.4.1          ffbase_0.12.7         
##   [7] bit64_0.9-7            knitr_1.33             aroma.light_3.14.0    
##  [10] R.utils_2.9.0          data.table_1.12.2      rpart_4.1-15          
##  [13] hwriter_1.3.2          RCurl_1.95-4.12        doParallel_1.0.14     
##  [16] generics_0.0.2         GenomicFeatures_1.36.4 preprocessCore_1.46.0 
##  [19] RSQLite_2.1.2          europepmc_0.3          bit_1.1-14            
##  [22] enrichplot_1.4.0       xml2_1.2.2             lubridate_1.7.4       
##  [25] assertthat_0.2.1       viridis_0.5.1          xfun_0.29             
##  [28] hms_0.5.2              jquerylib_0.1.4        evaluate_0.14         
##  [31] DEoptimR_1.0-8         fansi_0.4.0            progress_1.2.2        
##  [34] caTools_1.17.1.2       readxl_1.3.1           igraph_1.2.4.1        
##  [37] DBI_1.0.0              geneplotter_1.62.0     htmlwidgets_1.3       
##  [40] ellipsis_0.2.0         backports_1.1.4        bookdown_0.11         
##  [43] annotate_1.62.0        biomaRt_2.40.0         vctrs_0.2.0           
##  [46] withr_2.1.2            ggforce_0.2.2          triebeard_0.3.0       
##  [49] robustbase_0.93-5      checkmate_1.9.3        sna_2.4               
##  [52] prettyunits_1.0.2      cluster_2.0.8          DOSE_3.10.1           
##  [55] lazyeval_0.2.2         crayon_1.3.4           genefilter_1.66.0     
##  [58] pkgconfig_2.0.2        labeling_0.3           tweenr_1.0.1          
##  [61] nlme_3.1-139           nnet_7.3-12            rlang_0.4.2           
##  [64] lifecycle_0.1.0        modelr_0.1.5           cellranger_1.1.0      
##  [67] polyclip_1.10-0        graph_1.62.0           Matrix_1.2-17         
##  [70] urltools_1.7.3         base64enc_0.1-3        ggridges_0.5.1        
##  [73] viridisLite_0.3.0      bitops_1.0-6           R.oo_1.22.0           
##  [76] KernSmooth_2.23-15     blob_1.2.0             qvalue_2.16.0         
##  [79] robust_0.4-18          gridGraphics_0.4-1     ggsignif_0.5.0        
##  [82] scales_1.0.0           memoise_1.1.0          plyr_1.8.4            
##  [85] gdata_2.18.0           zlibbioc_1.30.0        compiler_3.6.0        
##  [88] intergraph_2.0-2       rrcov_1.4-7            cli_2.0.0             
##  [91] htmlTable_1.13.1       Formula_1.2-3          MASS_7.3-51.4         
##  [94] WGCNA_1.68             tidyselect_0.2.5       stringi_1.4.3         
##  [97] highr_0.8              yaml_2.2.0             GOSemSim_2.10.0       
## [100] locfit_1.5-9.1         latticeExtra_0.6-28    GeneOverlap_1.20.0    
## [103] grid_3.6.0             fastmatch_1.1-0        tools_3.6.0           
## [106] rstudioapi_0.10        foreach_1.4.4          foreign_0.8-71        
## [109] farver_1.1.0           digest_0.6.19          rvcheck_0.1.3         
## [112] pracma_2.2.5           ff_2.2-14              broom_0.5.3           
## [115] gRbase_1.8-3.4         httr_1.4.1             AnnotationDbi_1.46.0  
## [118] colorspace_1.4-1       rvest_0.3.5            XML_3.98-1.20         
## [121] splines_3.6.0          RBGL_1.60.0            ggplotify_0.0.3       
## [124] fit.models_0.5-14      xtable_1.8-4           jsonlite_1.6          
## [127] dynamicTreeCut_1.63-1  zeallot_0.1.0          R6_2.4.0              
## [130] Hmisc_4.2-0            pillar_1.4.3           htmltools_0.3.6       
## [133] glue_1.3.1             clusterProfiler_3.12.0 DT_0.7                
## [136] DESeq_1.36.0           codetools_0.2-16       pcaPP_1.9-73          
## [139] mvtnorm_1.0-11         lattice_0.20-38        network_1.15          
## [142] gtools_3.8.1           GO.db_3.8.2            survival_2.44-1.1     
## [145] rmarkdown_2.10.6       statnet.common_4.3.0   munsell_0.5.0         
## [148] DO.db_2.9              fastcluster_1.1.25     GenomeInfoDbData_1.2.1
## [151] iterators_1.0.10       impute_1.58.0          haven_2.2.0           
## [154] gtable_0.3.0